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We have developed a new Hartree-Fock-Bogoliubov (HFB) code which has been specificaUy de- 
signed to study ground state properties of nuclei near the neutron and proton drip lines. The unique 
feature of our code is that it takes into account the strong coupling to high-energy continuum states, 
up to an equivalent single-particle energy of 60 MeV. We solve the HFB equations for deformed, 
axially symmetric even-even nuclei in coordinate space on a 2-D lattice with Basis-Spline methods. 
" ".j", For the p-h channel, the Skyrme (SLy4) effective N-N interaction is utilized, and for the p-p and 

■ h-h channel we use a delta interaction. We present results for binding energies, deformations, nor- 

^-•^ ' mal densities and pairing densities, Fermi levels, and pairing gaps. In particular, we will discuss 

, neutron-rich isotopes of oxygen (^^O) and tin (^^"Sn). 
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I. INTRODUCTION 



One of the fundamental questions of nuclear structure physics is: what are the limits of nuclear stability? How many 
neutrons or protons can we add to a given nucleus before it becomes unstable against spontaneous neutron or proton 
emission? If one connects the isotopes with zero neutron separation energy, 5ri = 0, in the nuclear chart one obtains 
the neutron dripline. Similarly, the proton dripline is defined by the condition Sp ~ 0. Another limit to stability is 
the superheavy element region around Z — 124 — 126 and N = 184 which is formed by a delicate balance between 
0^ strong Coulomb repulsion and additional binding due to closed shells. The nuclear chart shows less than 300 stable 
, nuclear isotopes, and about 1700 additional isotopes have been synthesized and studied in accelerator experiments. 
Nuclei in between the proton and neutron driplines are unstable against /3-decay. Nuclei outside the driplines decay 
by spontaneous neutron emission or proton radioactivity. The neutron-rich side, in particular, exhibits thousands of 
nuclear isotopes still to be explored ('terra incognita'). 

Some of these exotic nuclei can be studied with existing first-generation Radioactive Ion Beam Facilities (e.g. HRIBF 
at Oak Ridge, NSCL at Michigan State University, GANIL in France, GSI in Germany, and RIKEN in Japan). Several 
^ countries are planning to construct new 'second generation' RIB facilities, in particular for the exploration of neutron 
rich isotopes. In the United States, the DOE/NSF Long Range Plan for Nuclear Physics, published in April 2002, 
ILj^ gives the highest priority for new construction to RIA (Rare Isotope Accelerator). RIA is a bold new concept in 
exotic beam facilities in that it combines both of the known rare isotope production techniques: the ISOL method 
(thick-target spallation) and high-energy projectile fragmentation. 
5^ Theories predict profound differences between the known isotopes near stability and the exotic nuclei at the driplines 

[0, 1^ : for n-rich nuclei, as the Fermi level approaches the particle continuum at = (see Fig. |l|) weakly bound states 
couple strongly to the continuum states giving rise to neutron halos and neutron skins. Theories also expect large 
pairing correlations and new collective modes (e.g. 'pigmy resonance'), a weakening of the spin-orbit force leading to 
a quenching of the shell gaps, and perhaps new magic numbers. 

Furthermore, RIA will allow us to address fundamental questions in nuclear astrophysics: more than half of all 
elements heavier than iron are thought to be produced in supernovae explosions by the rapid neutron capture process 
(r- process). The r-process path contains many exotic neutron- rich nuclei which can only be studied with RIA. Also, 
the predicted neutron skins would allow us to measure the properties of pure neutron matter which is of great interest 
for the study of neutron stars. 

These experimental developments as well as recent advances in computational physics have sparked renewed interest 
in nuclear structure theory. There are several types of approaches in nuclear structure theory [^: for the lightest 
nuclei, ab-initio calculations (Green's function Monte Carlo, no-core shell model) based on the bare N-N interaction 
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FIG. 1; One-body 'mean field' potential for protons and neutrons. Left: for a nucleus near the stability line. Right: for a 
nucleus near the neutron dripline. The Fermi levels for protons and neutrons are indicated by horizontal lines. 

are possible Medium-mass nuclei up to ^ ^ 60 may be treated in the large-scale shell model approach jij. For 
heavier nuclei one utilizes cither nonrelativistic ||^, ^ or relativistic 0, ||] mean field theories. 

II. HFB THEORY IN COORDINATE SPACE: DEFORMED NUCLEI ON A 2-D LATTICE 



An accurate treatment of the pairing interaction is essential for exotic nuclei 1^, ^. As we move away from the 
valley of stability, surprisingly little is known about the pairing force: for example, what is its density dependence? 
Neutron-rich nuclei are expected to be highly superfluid due to continuum excitation of neutron 'Cooper pairs'. These 
large pairing correlations near the driplines can no longer be described by a small residual interaction. It becomes 
necessary to treat the mean field and the pairing field in a single self-consistent theory, known as the Hartree-Fock- 
Bogoliubov theory (HFB). Near the neutron dripline, the Fermi energy approaches E — and the outermost nucleons 
are weakly bound (which implies a large spatial extent); they are strongly coupled to the particle continuum at i? > 0. 
These features represent major challenges for the numerical solution. 

Most HFB calculations to date are carried out in a truncated discrete harmonic ocillator basis, see e.g. |l^, 

III,!!!- 

This approach is quite appropriate for nuclei in the vicinity of the stability line. However, farther from stability, the 
continuum states become important and coordinate-based representations have numerous advantages: for example, 
the well-known 'French code' uses a truncated 3-D Hartree-Fock basis which consists of both localized states and 
discretized continuum states; however, in this approach one can only include continuum states up to about 5 MeV of 
excitation energy. For nuclei in the vicinity of the driplines, continuum states with an equivalent single-particle energy 
of up to 60 MeV must be taken into account. One-dimensional calculations for spherical nuclei have been carried out 
in coordinate space for many years |^ , but only recently has our Vanderbilt group succeeded in generalizing this 
approach to the important case of deformed axially symmetric nuclei (HFB on a 2-D lattice) . Wc utilize 

a novel computational technique, a Basis-Spline representation of wavefunctions and operators, which allows us to 
accurately describe high-energy continuum states in two space dimensions {z,r). 

The many-body Hamiltonian in occupation number representation has the form 



H^y<t\t\j> 



hi 
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The general linear transformation from particle operators c, to quasiparticle operators (3, (3^ take the form 



(2.2) 



The HFB approximate ground state of the many-body system is defined as a vacuum with respect to quasiparticles 

/3fc |$o > = . (2.3) 



3 



The HFB ground state energy including the constraint on the particle number N is given by 

E{n) ^<<i>a\H - XN\^o> (2.4) 

where the Lagrange multiplier A is the Fermi energy of the system. The equations of motion are derived from the 
variational principle 

S\n [E(n) - tr A(7^2 - 7^)] = , (2.5) 
where TZ represents the generalized density matrix. 

A. Quasiparticle wave functions and densities 

In practice, it is to convenient to transform the standard HFB equations into a coordinate space representation 
and solve the resulting differential equations on a lattice. For this purpose, one defines two types of quasiparticle 
wavefunctions (f>i and (f>2 

(f>UEc.,raq) = ^C/,„ (2a) </).(r-ag), MEo.,raq) = ^^/^(/..(mq). (2.6) 

i i 

The basis wavefunctions (j>i depend on the position vector r, the spin projection cr = ±i, and the isospin projection 
q (g = +i corresponds to protons and q = to neutrons). 

From these wavefunctions we obtain the following expressions for the normal density Pq(r) and the pairing density 
P«(r) 



Pq 



(r) = 5m'^2,a(rCT(?) 02,a(rCT'7)> Pg(i") = -5m'?^2.a(rcrg) <?!)*^„(rcrq) . (2.7) 



The quasiparticle energy Ea is denoted by index a for simplicity. In principle, the sums go over all the energy states, 
but in practice a cutoff is introduced (see later). The physical interpretation of pq has been discussed in the 
quantity [/5g(r) AV/2]'^ gives the probability to find a correlated pair of nucleons with opposite spin projection in the 
volume element AV. The kinetic energy density Tg(r) is found to be 

rqir)^J2J2\^ ■ (2-8) 



B. Binding energy functional 

In our calculations we utilize the Skyrme two-body effective N-N interaction 

v[l^ = to{l+ xoPa) 5{ri - ra) + i h (1 + xiP,) {5{ri - r2)k) + fc"(5(ri - ra)} 

+ t2 (1 + X2pa) k • '5(ri - ra) k + i i3 (1 + x^P„) 5{vi - V2) + i Wo (<ti + (72) ■ {k x 5(ri - r2)k}(2.9) 

The density-dependent term proportional to ^3 simulates the three-body N-N interaction v'^^\ 
The total binding energy of the nucleus 

Eq^^ — (^ol-ff I'I'o) — Ekin + Esky + Esky.LS + Ecoul + Epair + Ecm ■ (2-10) 

consists of a kinetic energy term, various contributions from the Skyrme effective N-N interaction (including a spin- 
orbit term). Coulomb and pairing energy, and a center-of mass correction due to the mean-field approximation 

2 
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(2.11) 



The Coulomb energy contains the direct term as well as an exchange term (in Slater approximation) 

Ecoui = y/^'^/ d'r'p,ir)j-l-j-^p,ir') - J d'r[p,{r)f' . (2.12) 



C. Pairing interaction 

In practice, one tends to use different effective N-N interactions for the p-h / h-p channels and for the p-p / h-h 
channels. Most pairing calculations utilize a local pairing interaction of the form 



Vpira,r' -a') = Vo (5(r - r') 6,^,, F(r) 



(2.13) 



This parameterization describes two primary pairing forces: a pure delta interaction (_F 1) that gives rise to volume 
pairing, and a density dependent delta interaction (DDDI) that gives rise to surface pairing. In the latter case, one 
uses the following phenomenological ansatz p3] for the factor F 



F(r) = 1- 



Po 



(2.14) 



where p(r) is the mass density. 

The pairing contribution to the nuclear binding energy is then 



F --^ 



I d'rY,p'(.r)F{r) 



(2.15) 



An important related quantity is the average pairing gap for protons and neutrons which can be calculated from the 
general expression given in |g, ^ 



<A,>= -i^y'dVp,(r)p,(r)F(r) 



(2.16) 



where Nq denotes the number of protons or neutrons. Note that the pairing gap is a positive quantity because Vq < 0. 

D. HFB equations and mean fields in coordinate space 

For certain types of effective interactions (e.g. Skyrme mean field and pairing delta-interactions) the particle 
Hamiltonian h and the pairing Hamiltonian h are diagonal in isospin space and local in position space, resulting in 
the following HFB equations with a 4x4 structure in spin space: 



(/i9 - A) h'i 
hi -(h^-X) 



with 



(2.17) 



(2.18) 



The HFB equations have a mathematical structure that is similar to the Dirac equation: the spectrum of quasiparticle 
energies E is unbounded from above and below. The spectrum is discrete for \E\ < —X and continuous for \E\ > —A. 
This is illustrated in Fig. |[ 
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FIG. 2: Computational challenge in solving the HFB equations in coordinate space: the quasiparticle energy spectrum is 
unbounded from above and below. 



As explained in it is forbidden to choose positive and negative quasiparticle energies at the same time, otherwise 
it is impossible to satisfy the anticommutation relations for Pa, Pa- For even-even nuclei it is customary to solve the 
HFB equations with a positive quasiparticle energy spectrum +Ea and consider all negative energy states as occupied 
in the HFB ground state. 

The HFB mean field Hamiltonian has the same structure as the binding energy functional 

K = -^'2;;^^ + + Uq.Lsir) + Ucouiir) ■ Sg.p . (2.19) 



The coordinate - dependent effective mass arises from the densities 

+ h p{r) - b[ p,{v) . (2.20) 



2m* {r) 2mq 



Detailed expressions for the Skyrme mean fields and the Coulomb term are given in reference jlTj. The DDDI 
interaction generates the following pairing mean field for the two isospin orientations g = ±^ 

h,{r) = i Vop,{r)F{r) . (2.21) 



III. 2-D REDUCTION FOR AXIALLY SYMMETRIC SYSTEMS 



For simplicity, we assume that the HFB quasiparticle Hamiltonian is invariant under rotations around the z-axis, 
i.e. [H,Rz] = 0. Due to the axial symmetry of the problem, it is advantageous to introduce cylindrical coordinates 
((/), r, z). It is possible to construct simultaneous eigenfunctions of the generalized Hamiltonian Ti, and the z-component 
of the angular momentum, jz with quantum numbers f2 = ±i,±|,±|,... corresponding to each nth energy state. 
The simultaneous quasiparticle eigenfunctions take the form 



'>Pn,n,qi4',r, z) 



We introduce the following useful notation 




(3.1) 



U 



(1,2) 



.(1,2) 



(r,z) = (/>i^,si,g('''^'T) , 

(r,z) = (t>n;n]g{r,z,l) . 



For axially symmetric systems, it is possible to eliminate the dependence on the angle 
problem in cylindrical coordinates [ pT| : 



(3.2) 
(3.3) 

resulting in the reduced 2-D 
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(3.4) 



with 




(1) ■ 



(r, z) 



n.fl.q 

(2)' 
n,fl.q 



(3.5) 



Here, quantities h', h' , U and L are all functions of {r, z) only. This is the main mathematical structure that 
we implement in computational calculations. For a given angular momentum projection quantum number Q, we 
solve the eigenvalue problem to obtain energy eigenvalues En.n^q and eigenvectors ipn,n,q for the corresponding HFB 
quasiparticle states. 

From the definitions of the normal density and pairing density we find the corresponding expressions in axial 
symmetry: 



Pqir, z) 




Pq 



X > I U%{r, z)U^,^^{r, z) + L^2^{r, z)i-;(r, z) 



(2) 



-(!)*, 



(3.6) 
(3.7) 



B,i>0 



IV. LATTICE REPRESENTATION OF SPINOR WAVEFUNCTIONS AND HAMILTONIAN 

We solve the HFB eigenvalue problem by direct diagonalization on a two-dimensional grid (r^jZ^), where a = 
1, ...,Nr and /3 = 1, ...,Nz. The four components of the spinor wavefunction are represented on the two-dimensional 
lattice by an expansion in basis-spline functions Bi{x) evaluated at the lattice support points. Further details about 
the basis-spline technique are given in Ref. ||l8|; |l9|. 

For the lattice representation of the Hamiltonian, we use a hybrid method |2^, |2^, ^ in which derivative operators 
are constructed using the Galerkin method; this amounts to a global error reduction. Local potentials are represented 
by the basis-spline collocation method (local error reduction). The lattice representation transforms the differential 
operator equation into a matrix form 

JV 

Y.n;i:^^E^^^ (^=l,...,iV) , (4.1) 

1^=1 

The calculations use as a starting point the result of a HF+BCS previous calculation, which makes HFB converge 
substantially faster. Since the problem is self-consistent we use an iterative method for the solution. At every iteration 
the full HFB Hamiltonian is diagonalized. Due to the axial symmetry in the intrinsic frame, the diagonalization is 
performed separately for each value of the angular momentum projection quantum number 17 and for the two isospin 
projections q = Typically 20-30 iterations are sufhcient for convergence at the level of one part in 10^ for the 
total binding energy. 

Note that in this lattice approach, the number of quasiparticle states is determined by the dimensionality of the 
discrete HFB Hamiltonian which is TV = (ANrN^)'^. 



V. NUMERICAL RESULTS 



In the following we discuss some numerical results obtained with our new HFB-2D code. First we present calculations 
for a light stable nucleus, fo^e. This nucleus was chosen because it has a large prolate g.s. quadrupole deformation. 
The calculation has been performed with the Skyrme SLy4 interaction in the p-h channel and a pure delta pairing 
interaction (strength Vb — — niMeVfm'^). For the pairing forces of zero range employed here, one needs to introduce 
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Total-Neutron-Potential Total-Proton-Potential 





FIG. 4: Neutron and proton density distribution for the strongly deformed nucleus Ne 



an energy cut-off. In all of our calculations, we use a cut-off energy Emax = 60 MeV in the equivalent s.p. energy 
spectrum^ the same value utilized by Dobaczewski et al. in his spherical 1-D calculations. 

Fig. ^ shows the mean field potential Uq{Y) for neutrons and protons. Both mean fields are fairly similar, except that 
there is an additional Coulomb contribution for the protons. In Fig. ^ we depict contour plots of the normal density 
for neutrons and protons. The large prolate quadrupole deformation is clearly visible in the density distributions. 
Fig. 1^ shows the pairing density for neutrons and protons. The square of the pairing density describes the probability 
of finding a correlated nucleon pair with opposite spin directions at position r Because this is a stable isotope, 
the pairing turns out to be relatively weak. The shell structure of this light nucleus causes significant differences in 
the neutron- vs. proton pairing. 
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FIG. 5: Neutron and proton pairing density distribution for Ne 



We now present results for the n-rich isotope ^^O which is close to the experimentally confirmed (2| dripline nucleus 
^'^O. Because this nucleus turns out to be spherical in our 2-D calculations, we can compare our results to the 1-D 
spherical code of Dobaczewski and to the 2-D oscillator basis expansion method of Stoitsov et al. . Table | 
shows several observables for this nucleus: the total binding energy, Fermi levels, pairing gaps and the r.m.s. radius. 
Overall, the results of the axially symmetric code of the present work agree with the other two in all observables. 



TABLE I: Comparison of several HFB calculations for ^^O. In all cases, the mean field is calculated with the SLY4 interaction, 
and the pairing force is a pure delta interaction (strength Vo = — 218.5Mel^/m'') corresponding to volume pairing. The axially 
symmetric calculations (2D) of this work used a box size R = 10/m with maximum jz = |. The spherical calculation of Ref. 
was made with R = 25 fm and a jmax = ^ . 



1-D 



25| 



2-D (TH0)[L3, 26 



2-D (this work) 



total binding energy (MeV) 
Fermi level, neutrons (MeV) 
Fermi level, protons (MeV) 
pairing gap, neutrons (MeV) 
pairing gap, protons (MeV) 
r.m.s. radius (fm) 



-164.60 
-5.26 

-18.88 
1.42 
0.0 
2.92 



-164.52 
-5.27 

-18.85 
1.41 
0.0 
2.92 



-164.64 
-5.29 
-18.08 
1.36 
0.0 
2.92 



We now turn our attention to the heavy neutron rich isotope ^^^Sn {Z = 50, N = 100) which has twice as many 
neutrons as protons. In some phenomenological models, it is already close to the 2-neutron dripline. Fig. ^ shows 
contour plots of the normal density for neutrons and protons. Because of the magic proton number and due to strong 
pairing, the nuclear shape turns out to be spherical in our 2-D HFB code. In Fig. |^ we depict a cut of these densities 
in radial direction which clearly exhibits the 'neutron skin' already seen in the 1-D calculations of Ref. The 
neutron skin represents a region of pure neutron matter which will allow us to the study the properties of neutrons 
stars in future laboratory experiments at RIA. 

In Table || we compare again our 2-D calculations with the 1-D radial results of Ref. ||] . In the 2-D calculations an 
approximately 3000 x 3000 matrix was diagonalized for each jz and isospin value, and for each major HFB iteration. 
The full calculation required about 20 HFB iterations. All observables agree quite well. The neutron Fermi level is 
less than IMeV in this case which demonstrates clearly the proximity of ^^^Sn to the two-neutron drip line. 
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Neutron-Density Proton-Density 




FIG. 6: Neutron and proton density distribution for foSn, a nucleus near the two-neutron dripline 



Sn (Z=50, N=100) 




r(fm) 

FIG. 7: Cut in radial direction of the neutron and proton density in ^loSn. The neutron skin is clearly visible. 

VI. SUMMARY AND CONCLUSIONS 

Our goal for the near future is to investigate several isotope chains, in particular deformed nuclei, and to calculate 
observables which are important for the physics near the drip lines, i.e. binding energies, neutron and proton separation 
energies, pairing gaps, particle densities and pairing densities, rms radii, and electric or magnetic moments. We plan 
to utilize a variety of Skyrme parameterizations for the mean field, and both volume and surface pairing forces. As 
more data from existing RIB facilities become available, it is likely that it will become necessary to develop new 
effective N-N interactions to describe these exotic nuclei. Furthermore, our 2-D HFB code results may be used as 
input into coordinate-space based QRPA calculations to investigate collective excited states (surface vibrations and 
giant resonances). 

We will also study alternative numerical techniques to speed up our 2-D HFB code, in particular damping methods 
which we have utilized successfully in solving the Dirac equation ]l9| . Provided that the damping method can be 
successfully implemented in 2-D, we will attempt to solve the 3-D HFB problem in Cartesian coordinates. This will 
certainly be very difficult, but it is worth trying: in the 1996 DOE/NSAC Long Range Plan, unrestricted HFB theory 
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FIG. 8: Neutron and proton pairing density distribution for foSn. Because of the magic proton number Z = 50 the proton 
pairing density is exactly zero. 



TABLE II: Comparison of 1-D and 2-D calculations for ^'""Sn. In both cases, the mean field is calculated with the SLY4 
interaction, and the pairing force is a pure delta interaction (strength Vb — —170MeVfm'^) corresponding to volume pairing. 
The 1-D calculations used a box size R = 30 with jmax of Calculations with our axially symmetric HFB 2-D code were 
made using a box size R = 20 fm with maximum 



Jmax 
13 
2 ■ 



1-D 



25] 



2-D (this work) 



total binding energy (MeV) 
Fermi level, neutrons (MeV) 
Fermi level, protons (MeV) 
pairing gap, neutrons (MeV) 
pairing gap, protons (MeV) 
r.m.s. radius (fm) 



-1129 
-0.96 
-17.54 
1.02 
0.00 
5.126 



-1130 
-0.94 
-17.34 
0.97 
0.00 
5.132 



on the lattice has been described as a computational Grand Challenge project in nuclear physics. 
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